package com.myApp.myaplicacion;

import java.util.Arrays;

import android.graphics.Bitmap;
import android.graphics.Bitmap.Config;

public class CannyEdgeDetector {

// statics

private final static float GAUSSIAN_CUT_OFF = 0.005f;
private final static float MAGNITUDE_SCALE = 10F;
private final static float MAGNITUDE_LIMIT = 100F;
private final static int MAGNITUDE_MAX = (int) (MAGNITUDE_SCALE * MAGNITUDE_LIMIT);

// fields

private int height;
private int width;
private int picsize;
private int[] data;
private int[] magnitude;
private Bitmap sourceImage;
private Bitmap edgesImage;

private float gaussianKernelRadius;
private float lowThreshold;
private float highThreshold;
private int gaussianKernelWidth;
private boolean contrastNormalized;

private float[] xConv;
private float[] yConv;
private float[] xGradient;
private float[] yGradient;

//add this to chase the stack, crash is at ~150 on Nexus One
private int mFollowStackDepth = 100; 

// constructors

/**
 * Constructs a new detector with default parameters.
 */

public CannyEdgeDetector() {
    lowThreshold = 2.5f;
    highThreshold = 3f;
    gaussianKernelRadius = 2f; //Por defecto 2f, valor alto menos ruido
    gaussianKernelWidth = 16; //Por defecto 16
    contrastNormalized = false;
}

// accessors

/**
 * The image that provides the luminance data used by this detector to
 * generate edges.
 * 
 * @return the source image, or null
 */

public Bitmap getSourceImage() {
    return sourceImage;
}

/**
 * Specifies the image that will provide the luminance data in which edges
 * will be detected. A source image must be set before the process method
 * is called.
 *  
 * @param image a source of luminance data
 */

public void setSourceImage(Bitmap image) {
    // Convert to RGB
    sourceImage = image.copy(Bitmap.Config.ARGB_8888, true);
    //sourceImage = image.copy(Bitmap.Config.ALPHA_8, true);
}

/**
 * Obtains an image containing the edges detected during the last call to
 * the process method. The buffered image is an opaque image of type
 * BufferedImage.TYPE_INT_ARGB in which edge pixels are white and all other
 * pixels are black.
 * 
 * @return an image containing the detected edges, or null if the process
 * method has not yet been called.
 */

public Bitmap getEdgesImage() {
    return edgesImage;
}

/**
 * Sets the edges image. Calling this method will not change the operation
 * of the edge detector in any way. It is intended to provide a means by
 * which the memory referenced by the detector object may be reduced.
 * 
 * @param edgesImage expected (though not required) to be null
 */

public void setEdgesImage(Bitmap edgesImage) {
    this.edgesImage = edgesImage;
}

/**
 * The low threshold for hysteresis. The default value is 2.5.
 * 
 * @return the low hysteresis threshold
 */

public float getLowThreshold() {
    return lowThreshold;
}

/**
 * Sets the low threshold for hysteresis. Suitable values for this parameter
 * must be determined experimentally for each application. It is nonsensical
 * (though not prohibited) for this value to exceed the high threshold value.
 * 
 * @param threshold a low hysteresis threshold
 */

public void setLowThreshold(float threshold) {
    if (threshold < 0) throw new IllegalArgumentException();
    lowThreshold = threshold;
}

/**
 * The high threshold for hysteresis. The default value is 7.5.
 * 
 * @return the high hysteresis threshold
 */

public float getHighThreshold() {
    return highThreshold;
}

/**
 * Sets the high threshold for hysteresis. Suitable values for this
 * parameter must be determined experimentally for each application. It is
 * nonsensical (though not prohibited) for this value to be less than the
 * low threshold value.
 * 
 * @param threshold a high hysteresis threshold
 */

public void setHighThreshold(float threshold) {
    if (threshold < 0) throw new IllegalArgumentException();
    highThreshold = threshold;
}

/**
 * The number of pixels across which the Gaussian kernel is applied.
 * The default value is 16.
 * 
 * @return the radius of the convolution operation in pixels
 */

public int getGaussianKernelWidth() {
    return gaussianKernelWidth;
}

/**
 * The number of pixels across which the Gaussian kernel is applied.
 * This implementation will reduce the radius if the contribution of pixel
 * values is deemed negligable, so this is actually a maximum radius.
 * 
 * @param gaussianKernelWidth a radius for the convolution operation in
 * pixels, at least 2.
 */

public void setGaussianKernelWidth(int gaussianKernelWidth) {
    if (gaussianKernelWidth < 2) throw new IllegalArgumentException();
    this.gaussianKernelWidth = gaussianKernelWidth;
}

/**
 * The radius of the Gaussian convolution kernel used to smooth the source
 * image prior to gradient calculation. The default value is 16.
 * 
 * @return the Gaussian kernel radius in pixels
 */

public float getGaussianKernelRadius() {
    return gaussianKernelRadius;
}

/**
 * Sets the radius of the Gaussian convolution kernel used to smooth the
 * source image prior to gradient calculation.
 * 
 * @return a Gaussian kernel radius in pixels, must exceed 0.1f.
 */

public void setGaussianKernelRadius(float gaussianKernelRadius) {
    if (gaussianKernelRadius < 0.1f) throw new IllegalArgumentException();
    this.gaussianKernelRadius = gaussianKernelRadius;
}

/**
 * Whether the luminance data extracted from the source image is normalized
 * by linearizing its histogram prior to edge extraction. The default value
 * is false.
 * 
 * @return whether the contrast is normalized
 */

public boolean isContrastNormalized() {
    return contrastNormalized;
}

/**
 * Sets whether the contrast is normalized
 * @param contrastNormalized true if the contrast should be normalized,
 * false otherwise
 */

public void setContrastNormalized(boolean contrastNormalized) {
    this.contrastNormalized = contrastNormalized;
}

// methods

public void process() {
    width = sourceImage.getWidth();
    height = sourceImage.getHeight();
    picsize = width * height;
    initArrays();
    readLuminance();
    if (contrastNormalized) normalizeContrast();
    computeGradients(gaussianKernelRadius, gaussianKernelWidth);
    int low = Math.round(lowThreshold * MAGNITUDE_SCALE);
    int high = Math.round( highThreshold * MAGNITUDE_SCALE);
    performHysteresis(low, high);
    thresholdEdges();
    writeEdges(data);
}

// private utility methods

private void initArrays() {
    if (data == null || picsize != data.length) {
        data = new int[picsize];
        magnitude = new int[picsize];

        xConv = new float[picsize];
        yConv = new float[picsize];
        xGradient = new float[picsize];
        yGradient = new float[picsize];
    }
}

//NOTE: The elements of the method below (specifically the technique for
//non-maximal suppression and the technique for gradient computation)
//are derived from an implementation posted in the following forum (with the
//clear intent of others using the code):
//  http://forum.java.sun.com/thread.jspa?threadID=546211&start=45&tstart=0
//My code effectively mimics the algorithm exhibited above.
//Since I don't know the providence of the code that was posted it is a
//possibility (though I think a very remote one) that this code violates
//someone's intellectual property rights. If this concerns you feel free to
//contact me for an alternative, though less efficient, implementation.

private void computeGradients(float kernelRadius, int kernelWidth) {

    //generate the gaussian convolution masks
    float kernel[] = new float[kernelWidth];
    float diffKernel[] = new float[kernelWidth];
    int kwidth;
    for (kwidth = 0; kwidth < kernelWidth; kwidth++) {
        float g1 = gaussian(kwidth, kernelRadius);
        if (g1 <= GAUSSIAN_CUT_OFF && kwidth >= 2) break;
        float g2 = gaussian(kwidth - 0.5f, kernelRadius);
        float g3 = gaussian(kwidth + 0.5f, kernelRadius);
        kernel[kwidth] = (g1 + g2 + g3) / 3f / (2f * (float) Math.PI * kernelRadius * kernelRadius);
        diffKernel[kwidth] = g3 - g2;
    }

    int initX = kwidth - 1;
    int maxX = width - (kwidth - 1);
    int initY = width * (kwidth - 1);
    int maxY = width * (height - (kwidth - 1));

    //perform convolution in x and y directions
    for (int x = initX; x < maxX; x++) {
        for (int y = initY; y < maxY; y += width) {
            int index = x + y;
            float sumX = data[index] * kernel[0];
            float sumY = sumX;
            int xOffset = 1;
            int yOffset = width;
            for(; xOffset < kwidth ;) {
                sumY += kernel[xOffset] * (data[index - yOffset] + data[index + yOffset]);
                sumX += kernel[xOffset] * (data[index - xOffset] + data[index + xOffset]);
                yOffset += width;
                xOffset++;
            }

            yConv[index] = sumY;
            xConv[index] = sumX;
        }

    }

    for (int x = initX; x < maxX; x++) {
        for (int y = initY; y < maxY; y += width) {
            float sum = 0f;
            int index = x + y;
            for (int i = 1; i < kwidth; i++)
                sum += diffKernel[i] * (yConv[index - i] - yConv[index + i]);

            xGradient[index] = sum;
        }

    }

    for (int x = kwidth; x < width - kwidth; x++) {
        for (int y = initY; y < maxY; y += width) {
            float sum = 0.0f;
            int index = x + y;
            int yOffset = width;
            for (int i = 1; i < kwidth; i++) {
                sum += diffKernel[i] * (xConv[index - yOffset] - xConv[index + yOffset]);
                yOffset += width;
            }

            yGradient[index] = sum;
        }

    }

    initX = kwidth;
    maxX = width - kwidth;
    initY = width * kwidth;
    maxY = width * (height - kwidth);
    for (int x = initX; x < maxX; x++) {
        for (int y = initY; y < maxY; y += width) {
            int index = x + y;
            int indexN = index - width;
            int indexS = index + width;
            int indexW = index - 1;
            int indexE = index + 1;
            int indexNW = indexN - 1;
            int indexNE = indexN + 1;
            int indexSW = indexS - 1;
            int indexSE = indexS + 1;

            float xGrad = xGradient[index];
            float yGrad = yGradient[index];
            float gradMag = hypot(xGrad, yGrad);

            //perform non-maximal supression
            float nMag = hypot(xGradient[indexN], yGradient[indexN]);
            float sMag = hypot(xGradient[indexS], yGradient[indexS]);
            float wMag = hypot(xGradient[indexW], yGradient[indexW]);
            float eMag = hypot(xGradient[indexE], yGradient[indexE]);
            float neMag = hypot(xGradient[indexNE], yGradient[indexNE]);
            float seMag = hypot(xGradient[indexSE], yGradient[indexSE]);
            float swMag = hypot(xGradient[indexSW], yGradient[indexSW]);
            float nwMag = hypot(xGradient[indexNW], yGradient[indexNW]);
            float tmp;
            /*
             * An explanation of what's happening here, for those who want
             * to understand the source: This performs the "non-maximal
             * supression" phase of the Canny edge detection in which we
             * need to compare the gradient magnitude to that in the
             * direction of the gradient; only if the value is a local
             * maximum do we consider the point as an edge candidate.
             * 
             * We need to break the comparison into a number of different
             * cases depending on the gradient direction so that the
             * appropriate values can be used. To avoid computing the
             * gradient direction, we use two simple comparisons: first we
             * check that the partial derivatives have the same sign (1)
             * and then we check which is larger (2). As a consequence, we
             * have reduced the problem to one of four identical cases that
             * each test the central gradient magnitude against the values at
             * two points with 'identical support'; what this means is that
             * the geometry required to accurately interpolate the magnitude
             * of gradient function at those points has an identical
             * geometry (upto right-angled-rotation/reflection).
             * 
             * When comparing the central gradient to the two interpolated
             * values, we avoid performing any divisions by multiplying both
             * sides of each inequality by the greater of the two partial
             * derivatives. The common comparand is stored in a temporary
             * variable (3) and reused in the mirror case (4).
             * 
             */
            if (xGrad * yGrad <= (float) 0 /*(1)*/
                ? Math.abs(xGrad) >= Math.abs(yGrad) /*(2)*/
                    ? (tmp = Math.abs(xGrad * gradMag)) >= Math.abs(yGrad * neMag - (xGrad + yGrad) * eMag) /*(3)*/
                        && tmp > Math.abs(yGrad * swMag - (xGrad + yGrad) * wMag) /*(4)*/
                    : (tmp = Math.abs(yGrad * gradMag)) >= Math.abs(xGrad * neMag - (yGrad + xGrad) * nMag) /*(3)*/
                        && tmp > Math.abs(xGrad * swMag - (yGrad + xGrad) * sMag) /*(4)*/
                : Math.abs(xGrad) >= Math.abs(yGrad) /*(2)*/
                    ? (tmp = Math.abs(xGrad * gradMag)) >= Math.abs(yGrad * seMag + (xGrad - yGrad) * eMag) /*(3)*/
                        && tmp > Math.abs(yGrad * nwMag + (xGrad - yGrad) * wMag) /*(4)*/
                    : (tmp = Math.abs(yGrad * gradMag)) >= Math.abs(xGrad * seMag + (yGrad - xGrad) * sMag) /*(3)*/
                        && tmp > Math.abs(xGrad * nwMag + (yGrad - xGrad) * nMag) /*(4)*/
                ) {
                magnitude[index] = gradMag >= MAGNITUDE_LIMIT ? MAGNITUDE_MAX : (int) (MAGNITUDE_SCALE * gradMag);
                //NOTE: The orientation of the edge is not employed by this
                //implementation. It is a simple matter to compute it at
                //this point as: Math.atan2(yGrad, xGrad);
            } else {
                magnitude[index] = 0;
            }
        }
    }
}

//NOTE: It is quite feasible to replace the implementation of this method
//with one which only loosely approximates the hypot function. I've tested
//simple approximations such as Math.abs(x) + Math.abs(y) and they work fine.
private float hypot(float x, float y) {
    return (float) Math.hypot(x, y);
}

private float gaussian(float x, float sigma) {
    return (float) Math.exp(-(x * x) / (2f * sigma * sigma));
}

private void performHysteresis(int low, int high) {
    //NOTE: this implementation reuses the data array to store both
    //luminance data from the image, and edge intensity from the processing.
    //This is done for memory efficiency, other implementations may wish
    //to separate these functions.
    Arrays.fill(data, 0);

    int offset = 0;
    for (int x = 0; x < width; x++) {
        for (int y = 0; y < height; y++) {
            if (data[offset] == 0 && magnitude[offset] >= high) {
                follow(x, y, offset, low, mFollowStackDepth);
            }
            offset++;
        }
    }
}

private void follow(int x1, int y1, int i1, int threshold, int depth)
{
       if( depth > mFollowStackDepth)  // don't run out of stack!
           return;
       int x0 = x1 == 0 ? x1 : x1 - 1;
       int x2 = x1 == width - 1 ? x1 : x1 + 1;
       int y0 = y1 == 0 ? y1 : y1 - 1;
       int y2 = y1 == height -1 ? y1 : y1 + 1;

       data[i1] = magnitude[i1];
       for (int x = x0; x <= x2; x++) {
           for (int y = y0; y <= y2; y++) {
               int i2 = x + y * width;
               if ((y != y1 || x != x1)
                   && data[i2] == 0 
                   && magnitude[i2] >= threshold) {
                   follow(x, y, i2, threshold, depth+1);
                   return;
               }
           }
       }
   }

private void thresholdEdges() {
    for (int i = 0; i < picsize; i++) {
        data[i] = data[i] > 0 ? -1 : 0xff000000;
    }
}

private int luminance(float r, float g, float b) {
    return Math.round(0.299f * r + 0.587f * g + 0.114f * b);
}

private void readLuminance() {
    //int type = sourceImage.getType();
    //if (type == BufferedImage.TYPE_INT_RGB || type == BufferedImage.TYPE_INT_ARGB) {
    // We short-circuit this because we manually set the image as ARGB earlier
    if (true)
    {
        //int[] pixels = (int[]) sourceImage.getData().getDataElements(0, 0, width, height, null);
        int[] pixels = new int[picsize];
        sourceImage.getPixels(pixels, 0, width, 0, 0, width, height);
        for (int i = 0; i < picsize; i++) {
            int p = pixels[i];
            int r = (p & 0xff0000) >> 16;
            int g = (p & 0xff00) >> 8;
            int b = p & 0xff;
            data[i] = luminance(r, g, b);
        }
    }
    /*
    elseif (type == BufferedImage.TYPE_BYTE_GRAY) {
        //byte[] pixels = (byte[]) sourceImage.getData().getDataElements(0, 0, width, height, null);
        int[] pixels = new int[picsize];
        sourceImage.getPixels(pixels, 0, width, 0, 0, width, height);
        for (int i = 0; i < picsize; i++) {
            data[i] = (pixels[i] & 0xff);
        }
    }
    /*else if (type == BufferedImage.TYPE_USHORT_GRAY) {
        short[] pixels = (short[]) sourceImage.getData().getDataElements(0, 0, width, height, null);
        for (int i = 0; i < picsize; i++) {
            data[i] = (pixels[i] & 0xffff) / 256;
        }
    } else if (type == BufferedImage.TYPE_3BYTE_BGR) {
        byte[] pixels = (byte[]) sourceImage.getData().getDataElements(0, 0, width, height, null);
        int offset = 0;
        for (int i = 0; i < picsize; i++) {
            int b = pixels[offset++] & 0xff;
            int g = pixels[offset++] & 0xff;
            int r = pixels[offset++] & 0xff;
            data[i] = luminance(r, g, b);
        }
    } else {
        throw new IllegalArgumentException("Unsupported image type: " + type);
    }
    */
}

private void normalizeContrast() {
    int[] histogram = new int[256];
    for (int i = 0; i < data.length; i++) {
        histogram[data[i]]++;
    }
    int[] remap = new int[256];
    int sum = 0;
    int j = 0;
    for (int i = 0; i < histogram.length; i++) {
        sum += histogram[i];
        int target = sum*255/picsize;
        for (int k = j+1; k <=target; k++) {
            remap[k] = i;
        }
        j = target;
    }

    for (int i = 0; i < data.length; i++) {
        data[i] = remap[data[i]];
    }
}

private void writeEdges(int pixels[]) {
    //NOTE: There is currently no mechanism for obtaining the edge data
    //in any other format other than an INT_ARGB type BufferedImage.
    //This may be easily remedied by providing alternative accessors.
    if (edgesImage == null) {
        //edgesImage = new BufferedImage(width, height, BufferedImage.TYPE_INT_ARGB);
        edgesImage = Bitmap.createBitmap(width, height, Config.ARGB_8888);
    }
    //edgesImage.getWritableTile(0, 0).setDataElements(0, 0, width, height, pixels);
    edgesImage.setPixels(pixels, 0, width, 0, 0, width, height);
}
}